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ABSTRACT 

Context, (abridged) The physics of the pulsar magnetosphere near the neutron star surface re- 
mains poorly constrained by observations. Indeed, little is known about their emission mecha- 
nism, from radio to high-energy X-ray and gamma-rays. Nevertheless it is believed that large 
vacuum gaps exist in this magnetosphere, and a non-neutral plasma partially fills the neutron star 
surroundings to form an electrosphere in differential rotation. 

Aims. According to several of our previous works, the equatorial disk in this electrosphere is 
diocotron and magnetron unstable, at least in the linear regime. To better assess the long term 
evolution of these instabilities, we study the behavior of the non-neutral plasma with help on 
particle simulations. 

Methods. We designed a 2D electrostatic PIC code in cylindrical coordinates, solving Poisson 
equation for the electric potential. In the diocotron regime, the equation of motion for particles 
obeys the electric drift approximation. As in the linear study, the plasma is confined between two 
conducting walls. Moreover, in order to simulate a pair cascade in the gaps, we add a source term 
feeding the plasma with charged particles having the same sign as those already present in the 
electrosphere. 

Results. First we checked our code by looking for the linear development of the diocotron in- 
stability in the same regime as the one used in our previous work, for a plasma annulus and for 
a typical electrosphere with differential rotation. To very good accuracy, we retrieve the same 
growth rates giving confidence for the correctness of our PIC code. Next, we consider the long 
term non-linear evolution of the diocotron instability. We found that particles tend to attract to- 
gether to form small vortex of high charge density rotating around the axis of the cylinder with 
only little radial excursion of the particles. This grouping of particles generates new low den- 
sity or even vacuum gaps in the plasma column. Finally, in more general initial configurations, 
we show that particle injection into the plasma can drastically increase the diffusion of particles 
across the magnetic field lines. Also, it has to be noted that the newly formed vacuum gaps cannot 
be replenished by simply invoking the diocotron instability. 
Conclusions. 
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1. INTRODUCTION 



The detailed structure of charge distribution and electric-cuiTent circulation in the closed magne- 
tosphere of a pulsar remains puzzling. Although it is often assumed that the plasma fills the space 
entirely and corotates with the neutron star, it is on the contrary very likely that it only partly 
fills it, leaving large vacuum gaps between plasma- filled regions. The existence of such gaps in 
aligned rotators has been very clearly established by 



Krause-Polstorff & Michel! (llQSsXl) . Since 



then, a number of diff'erent n umerical app r oaches to the p roblem have confirmed their conclu 



sions, including some wo r kbvlRvlovl(ll989h.lShibatal(ll989l).IZachariadesl (ll993h.lNeukircl] (11993 ), 



Thielheim & Wolfstellej (1 1994 . 



Spitkovskv & AronsI (120021) . and ourselves (IPetri et al. 



2002b). 



This conclusion about the existence of vacuum gaps has been reached from a self-consistent so- 
lution of the Maxwell equations in the case of the aligned rotator. Moreover, 



Smith et al 



(120011) 

have shown by numerical modelling that an initially filled magnetosphere like the Goldreich- Julian 
model evolves by opening up large gaps and stabili zes to the pa rtially filled and pa rtially void so- 

J2002bl) . The status of 



Petri et al 



lution found bv lKrause-Polstorff" & Michell (Il985al) and also by 
models of the pulsar magnetospheres, or electrospheres, has recently been critically reviewed by 



Michell ( l2005b . A solution with vacuum gaps has the peculiar property that those parts of the mag- 



netosphere that are separated from the star's surface by a vacuum region are not corotating and so 
suffer differential rotation, an essential ingredient that will lead to non-neutral plasma instabilities 
in the closed magnetosphere. 

This kind of non-neutral plasma insta bility is well known to plasma physicists dOneill ( Il980h : 



Davidson 



(ll990l) : IO'Neil & SmithI (Il992h ). Their good confinement properties makes them a valu- 
able tool for studying plasmas in laboratory, by using for instance Penning traps. In the magneto- 
sphere of a pulsar, far from the light cylinder and close to the neutron star surface, the instability 
reduces to its non-relativistic and electrostatic form, the diocotron instability. The linear develop - 



ment of this instability for a d ifferentially ro tating charged disk was studied by 



in the thin disk limit, and by iPetril (l2007albh in the thick disk limit. It both cases, the instability 



Petri et al. 



(l2002ah . 



proceeds at a growth rate co mparable to the s tar's rotation rate. The non-linear development of 
this instability was studied by 



Petri et al 



(120031) . in the framework of an infinitely thin disk model. 
They have shown that the instability causes a cross- field transport of these charge s in the equatorial 
disk, evolving into a net out-flowing flux of charges. ISpitkovskv & AronsI ( l2002h have numerically 
studied the problem and concluded that this charge transport tends to fill the gaps with plasma. 
The appeara nce of a cross-field elect ric current as a result of the diocotron instability has been 
observed bv lPasquini & FaiansI (|2002|) in laboratory experiments in which charged particles were 
continuously injected in the plasma column trapped in a Malmberg-Penning configuration. 

Numerical simulations of the non-linear development of th e diocotron instability have been in- 



vestigated by different authors. For instance, 



Petri et al 



(I2OO3I) used a numerical scheme similar to 



those used for MHD codes. This has the drawback to hardly handle vacuum gaps which eventu 
afly are created in the electrosphere. To better deal with these gaps, we decided to design a parti 
cle simulation code that does not suffer from the prese nce of vacuum. PIC me t hods have already 
been s uccessful! to investigate the diocotron instability dNeu & MoralesI (ll995 K 



(I2OO3I) ). Problems with injection have also been considered, see for instance 



Yatsuvanagi et al 



Bettega et al 



(2007). 
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In this paper, we extend the work initiated bv lPetril ( l2007al) and look for the non-linear evolution 
of the diocotron instability by performing 2D electrostatic PIC simulations. The question on the 
influence of pair creation will also be addressed by permitting injection of charged particles into 
the system. The paper is organ i zed as f ollows. In Sect.|2] we briefly recall the initial setup of the 
plasma column as done in lPetril (l2007bl) . In Sect. [3] we show the results of our 2D PIC simulations. 
First we checked the code by computing some linear growth rate for special cases l ike a constan t 
density plasma column for which we know exact analytical solutions, see for instance lPetril (l2007al) . 



S econd , we che ck that we retrieve the same linear growth rates for the electrosphere as those found 
in lPetril ( l2007bl) . Third, we let the system evolve on long time scale and look for significant particle 



difiiision. Finally, in a last set of runs, we injected some charged particles into the system to study 
the behavior of the instability in the presence of a source of charges. Sect. |4| The extreme case 
starting with vacuum will also be presented. The conclusions and the possible generalizations are 
presented in Sect.|5] 



2. THE MODEL 



Let us first remind the plasma configuration as described in lPetril ( l2007bl) . We study the motion of 
a non-neutral plasma column of infinite axial extend along the z-axis. We adopt cylindrical coor- 
dinates denoted by (r,(p,z) and define the corresponding orthonormal basis vectors by (Cr,e^,ez). 
In the initial state corresponding to an equilibrium, the plasma is only rotating in the azimuthal 
direction, with no motion in the radial direction along r or vertical direction along z. 



2. 1 . Plasma and field setup 

Unlike other magnetospheric models, our pulsar electrospheric plasma is assumed to be completely 
charge-separated. So we consider a single-species non-neutral fluid consisting of particles with 
mass nie and charge q trapped between two cylindrically conducting walls located at the radii Wi 
and W2 > W]. The plasma column itself is confined between the radii Ri > Wi and R2 < W2, 
with R\ < R2. This allows us to take into account vacuum regions between the plasma and the 
conducting walls. Such geometric configuration will also be useful to investigate particle diffusion 
across magnetic field lines into this vacuum and will clearly demonstrate the efficiency of the 
diocotron instability to fill gaps with plasma as shown in some examples in the last section. 

In order to simulate the presence of a charged neutron star generating a radial electric 
quadrupole field coming from the rotating magnetic dipole, the inner wall at Wi can carry a 
charge Q per unit length such that its electric field is simply given by Maxwell-Gauss theorem, 
namely, (we use MKSA units) 

EJr)^-^e,. (1) 
ZTTSor 

In the equilibrium configuration, the particle number density is n(r) and the charge density 
is Pe{r) - qn{r). In contrast to earlier studies, the external magnetic field, along the z-axis, is 
not necessarily uniform but can decrease with radius 

B^BAr)e^. (2) 

The electric field is made of two parts, the first one arising from the plasma column itself, 
and the second one from the inner conducting wall E^, Eq. ([Tji. We assume that the electric field 
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induced by the plasma vanishes at r = Wi, i.e. Ep(Wi) = 0. Applying Maxwell-Gauss equation, 
we solve for the equilibrium electric field and get 

Ep(r)^— [ pjr') r' dr' e,. (3) 

The total electric field, directed along the radial direction is therefore the sum 

E^Ep+E„ = Ere,. (4) 

In the non-relativistic diocotron regime, we do not solve the full set of Maxwell equations but 
only the electrostatic part. This means that the magnetic field perturbation induced by the plasma 
flow is neglected. Thus, the externally imposed magnetic field, Eq. (|2]), remains constant during the 
simulations. This assumption is justified when the plasma is confined well inside the light-cylinder 
and when its rest mass energy density remains negligible compared to the magnetic field energy 
density, in other words 

2 

nnieC <K (5) 

2.2. Equation of motion 

The motion of the plasma is governed by the conservation of charge, which in the non-neutral 
plasma case is equivalent to conservation of the number of particles, the Maxwell-Poisson equation, 
and the electric drift approximation, respectively 

^+div(peV) = (6) 

ot 

+ ^ = (7) 

so 

E aB 

v-^ (8) 
E = -V0. (9) 

(p is the usual electric potential associated with the total electric field, Eq. (HJl (plasma distribution 
+ externally appUed electric field). We introduced the usual notation, pe for the electric charge den- 
sity, V for the velocity, {E, B) for the electromagnetic field. The set of Eqs. (|6])-(|9| fully describes 
the non-linear time evolution of the cold plasma in the non-relativistic diocotron regime. Note that 
particle inertia do not enter into the above mentioned expressions. The gyromotion is averaged 
over the fastest timescale, leading to the electric drift approximation, as long as this drift speed 
remains smaller than the speed of light, in other words, as long as £ < c B. In the above mentioned 
equations (|6])-(|9]l, there is no reference to the speed of light because we are in the non-relativistic 
regime. Formally it corresponds to the limit c — > -i-oo, electromagnetic waves propagate instan- 
taneously. Therefore, at this stage, it is impossible to check whether the electric field intensity E 
remains smaller than cB or: not. The hypothesis E < c B is implicitly always true in our model 
because of the electric drift approximation! Indeed v < c in Eq. ([8]) implies E < cB. It is an im- 
plicit assumption that we presume to be valid during the whole simulations. To include relativistic 
effects, we would have to solve the full set of Maxwell equations as was already done to inves tigate 
the rela tivistic stabilization of the diocotron instability in the linear regime, see for instance 



Petri 



( l2007ah . This would introduce a timescale related to the speed of Ught c. Performing 2D relativistic 
PIC simulations of the diocotron instability is the purpose of ongoing work and will be exposed in 
a forthcoming paper 
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2.3. Code description 



We designed an algorithm usi ng standard techniques for particle in cell (PIC) codes as detailed 



2005 



Hocknev R.W. 



1988h . We use expUcit 



for instance in two textbooks dBirdsall & Langdoni 
schemes to advance in time both the particles and the electric field. A description of our own code 
on how to evolve field and particles is briefly exposed in the following subsections. 



2.3.1 . Integration of the equation of motion 

We assume that the inner and outer conductors act mechanically as solid walls in such a way that 
particles are reflected at these boundaries. The equation of motion for particles in the electric drift 
approximation takes a very simple form allowing to solve immediately for the velocity vector, see 
Eq. ([8]l. From a computational point of view, this equation of motion, Eq. ([8]), is integrated with 
a leapfrog scheme, which is second order in time. More precisely, denoting t" = n Af the time 
at the n-th iteration where Af is the time step, positions are computed at integer times t" whereas 
velocities are computed at half-integer times f"+i/2 = (^n + 1/2) At for each particle labelled by the 
index / such that 

= vr"' (11) 



Af 

The electromagnetic field at half-integer time is evaluate from the position of the particles at the 
same time, positions updated according to the speed known only at preceding times f" '^^ like 

^„+l/2^ „ Af 

2 ' 

Then, the electromagnetic field at half-integer time following from these positions at the location 
of particle number / are 

£«+l/2 ^ £(^«+l/2) (13) 
g«+l/2 ^ B(^«+l/2) (14) 

We use a Lagrangian description which means that no grid is associated with the position of the 
particle. Because the initial velocity and position are given for f = by v'' and rj*, we initiahze the 
speed of the particles by first doing half a time step forwards in order to find 

rl'^-r^^^v^ (15) 
1/2 Ef^B]'' 

(b;/2) 

We then evolve positions and velocities according to Eqs. ( fTOl i. (fTTT i. Note that in all the simulations 
the electric and magnetic fields do not explicitly depend on time. 



2.3.2. Field solver 

Poisson equation is solved on an Eulerian grid by Fourier transform in the azimuthal direction ip 
and a finite difference scheme for the radial dependence. Note that in Eq. (fTOt both the electric and 
the magnetic fields are evaluated at a half-integer time f"+'^^. Because the electric field is derived 
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from a potential 0, we need only to consider a scalar Poisson equation. The potential (p as well as 
the charge density are developed into real cos / sin Fourier series as 



(f>(r, (p,t) ^ ^ fjr, t) cos(m tp) + (f>J„(r, t) sin(TO ip) (17) 

m=0 

pAr, ip,t) = ^ p4(r, f) cos(ot ip) + p,'^^{r, t) sin(m ip) (18) 



m=0 



A^^ is the number of discretization points in the azimuthal direction. Thus Poisson equation trans- 
forms into a set of ordinary differential equations for the unknown functions (j)m '\r) as follows 



1 d I ^^'VM „.2 „ cIs 



r dr 



J./ 



r ■ 



dr 



Because of the linearity of the Poisson equation, the azimuthal modes decouple into a set of inde- 
pendent ordinary differential equations for each natural integer m. Using a classical finite difference 
discretization in radius, we derive a set of matrices for each couple of unknown functions (0J„, <^^) 
to be solved by standard linear algebra techniques. However, the structure of the discretization 
matric es is tridiagonal allowing to use numerical algorithm that are more efficient, see for in- 



stance 



Press W. 



( 120021) . Moreover, as outer boundary condition, we insure a vanishing electric 
potential at the outer wall which is expressed as '(^2) = whereas at the inner boundary, we 
enforce either a vanishing derivative d(p^Ji\Wi)/dr - corresponding to nullifying the normal com- 
ponent of the electric field for the vacuum or plasma annuli cases, or the corotation electric field 
for the electrosphere, in order to adjust smoothly to the revolution of the compact object (assuming 
ideal MHD). 

Charges are deposited on the Eulerian grid (r, p) in order to deduce the source terms p^m ' 
for Poisson equation. The electromagnetic field is evaluated on particle positions by 2D linear 
interpolation between points in this Eulerian grid (r, ip). 



2.3.3. Particle injection 

Finally, we simulate the pair creation mechanism in the magnetosphere by injecting charged par- 
ticles somewhere inside the system. The deposition is uniform in the azimuthal direction and pre- 
scribed to be at a given radius and rate, both fixed by the user. Typically we took an injection 
radius 7?inj -2.0Wi. Formally, the shape function can be though as a delta-Dirac distribution in the 
radial direction, let us say S{r, tp) - 2n rT 6{r - /?inj). T is the particle flux. Nevertheless, this is 
a user defined source function and can be taken to behave very differently, that can also evolve in 
space and time. There is no particular restriction to its shape. As a starting point, we took the most 
simple dependence, constant in time and ip and sharply concentrated in radius. To be more realistic, 
pair creation should be incorporate by taking into account different processes Uke photon-photon 
disintegration or photon-(strong) magnetic field interaction (which is formally the same process as 
the former, as the magnetic field can be though as a photon field from a quantum mechanical point 
of view). Monte Carlo simulations would be a valuable tool for such a study, but this will not be 
included in this work. 
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3. RESULTS WITHOUT PARTICLE INJECTION 



7 



First, we checked our algorithm in different configurations by picking out the linear regime from 
our runs and compare the growth rates with those predicted by the linear analysis. 

In order to keep the particle noise as weak as possible to really see the linear stage evolving 
on a few order of magnitudes, we need to start with an initial state with minimal energy. This is 
achieve d by using a so-called "quie t start" such that space is regularly filled with particles, see for 



mstance 



Birdsall & Langdonl(l2005l) . 



The resolution of the spatial grid on which the electric potential is solved and charge is de- 
posited is A^, X Ntp - 200 X 256, if not otherwise specified. On average, at the beginning of the 
runs, there are 25 particles per cell giving a total number of 1,280,000. Time is given in units of 
l/a>B where cjb = qBlnig is the cyclotron frequency. We also tried different resolutions and find 
no qualitative discrepancy. So we adopted the canonical values aforementioned. 

Note that for all the runs shown below, the charge Q carried by the inner conductor is set to 
zero, Q-Q. The purpose of this work is not to make a parametric study of the diocotron instability 
behavior but rather to point out its main non-linear characteristics for a closed (without particle 
injection or loss) and an open (with particle injection or losses) system. Moreover, in this and the 
subsequent sections, the magnetic field is uniform and constant in the whole plasma column. Let 
us now discuss the main results. 



3.1. Plasma annulus 



For completeness, we briefly recall the main characteristics of the configuration. The magnetic 
field is constant and uniform outside the plasma annulus, namely = Bq, for r > Rj. We solve 
Maxwell-Poisson equation in the space between the inner and the outer wall, R\ < r < R2. The 
particle number density and charge density are constant in the whole plasma column such that 



, <r<Ri 

po = noq - const , Ri < r < R2 
, R2<r<W2 



(20) 



To initiate the instability we disturb the axisymmetric configuration by overlapping a tiny density 
perturbation with relative magnitude (i.e. compared to the background density) of the order h - 
10"^. To pick out a given mode, the perturbation is modulated with the appropriate azimuthal 
pattern m. To be more precise, the initial axisymmetric equihbrium position of each particle (ro, (fo) 
is perturbed according to 

V'init ^fa + h cos(m (po) (21) 

In this cylindrical geometry, an ex act analytical sol ution of the dispersion relation has been 
found for the diocotro n instability, s ee 
as was already done in 



Petri 



DavidsonI (119901) . We use these results to check our code. 



(l2007ah . Some exact analytical eigenvalues are listed in Table [T] for 
different radial plasma extensions and azimuthal modes. In order to quantify the growth rate of the 
diocotron instability in our simulations, we estimate the total electrostatic energy in the system by 
computing the integral 



WAt) 



Jw, Jo 



■ rdrdip 



(22) 
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Mode m 




d2 


0) 




2 


0.4 


0.5 


-3.772e-01 +7.176e-02 


i 


3 


0.4 


0.5 


-5.456e-01 + 2.267e-01 


i 


4 


0.4 


0.5 


-7.216e-01 + 2.988e-01 


i 


5 


0.7 


0.9 


-1.147e+00 + 5.787e-02 


i 


7 


0.6 


0.7 


-9.315e-01 + 3.307e-01 


i 



Table 1. Eigenvalues O) for the constant density plasma annulus, for different modes m and dif- 
fe rent aspect ratios , di = R1/W2, and d2 - R2/W2, obtained from the analytical expression found 



Davidson 



(1990. 



To clearly point out the linear growing phase, we calculate the increase in electrostatic energy by 
seeking for the difference between the actual and the initial electrostatic energy as 

AWeit) = We(t) - WJO) (23) 

Note that because AWdt) is a quadratic expression with respect to the perturbed quantities, the 
linear growth rate for the electrostatic energy increase will be twice the value given in Table [T] 

Two examples of run are now discussed. For the first run, the plasma annulus is located between 
/?i = 4 and R2 - 5 with a perturbation m = 3, call it Al while for the second run we took Ri - 6, 
R2—I with a perturbation m - 7, call it A2. Let us have a look on the evolution of the electrostatic 
energy given by Eq. ( |23] ) and depicted in Fig.[T]on a logarithmic scale. 

For the first run, at small times, let say f < 100 and after a while, the linear regime sets in 
and the measured growth rate is close to the one predicted by Unear analysis (red straight line), the 
perturbed pattern grows at a speed predicted by the linear theory. The plasma column stays roughly 
confined between /?i and R2, no strong density rearrangement is perceptible. Note also that no other 
mode is excited, no cascade due to non-linear effects is observed because to weak. 

In some runs, as the first aforementioned, it can happen that the initial electrostatic fluctuations 
decrease sUghtly before a significant increase of many orders of magnitude, after a while, inducing 
a AWeif) < at the beginning and thus making the logarithm not defined at this early stage. In 
FigH it is doubtful to give a physical interpretation to this initial behaviour because it corresponds 
to numerical round off error, close to the machine precision (we work with double precision num- 
bers, i.e. 16 significant digits). Moreover, the gaps can also be are probably due to the particle noise 
inherent of the PIC method, proportional to the inverse of the square root of the number of parti- 
cles A^, (i.e. oc N^^^^). Tiny fluctuations in the positions, especially particles receding to the inner 
wall, tend to decrease the total electrostatic energy, at least at the beginning of the simulations. 

In a second step, non-linear effects come into play and the initial plasma distribution is de- 
stroyed. The azimuthal pattern with the typical m = 3 mode becomes clearly visible, left part of 
Fig.|2] Particles "merge" together to form large vortices and leave behind them vacuum gaps which 
will not be replenished. For some time, these vortices rotate around the axis of the cylinder keeping 
their m = 3 structure. In a last step, almost all particles hit the inner wall W\, this happens when 
t > 400. The total electrostatic energy saturates, reaching a plateau at the latest time. When the 
non-linear effects set in, we are interested in the diffusion of particle across magnetic field lines 
and on its consequences on a longer timescale. To get a better idea of this diffusion process, we 
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Electrc^static energy Electrc^static energy 




100 200 JOO 400 500 600 100 200 JOO 400 500 600 

t t 



Fig. 1. Evolution of the total electrostatic increase in energy AWe(t) in the constant density plasma 
column obtained from our PIC simulation (black dots) and compared with the linear growth rate 
of the diocotron instability obtained from Tab. [1] (red line), on the left for m - 3, on the right for 
m-1. The linear relation holds for more than ten orders of magnitude. 



t = 12 5 . 664 t = 81 . 6814 




Fig. 2. Snapshot of the charge density in the plasma column showing the m - 3 pattern (on the 
left) and the m - 7 pattern (on the right). The chosen time corresponds to the transition between 
the linear phase and the beginning of the non-linear regime, associated with the total electrostatic 
energy curves discussed in Fig.[T] 



quantify it by computing the azimuthally integrated charge density as 



p2n 

Nir,t)= p,rd<p (24) 
Jo 

It can be interpreted as an equivalent radial charge density. Indeed 

N(r,t)dr (25) 

gives the total charge contained in the plasma column at time t. Some examples related to the 
simulations shown in Fig.[T]are given in Fig.|3] We immediately see that for t < 100 (linear stage), 
the boundary of the plasma column remains unchanged as discussed before. In a second stage, non- 
linear effects become important and strongly disturb the boundaries of the plasma annulus. Some 
very coherent patterns emerge while continuing to rotate. During this period roughly 100 < f < 400, 
the barycenter of the charges moved closer to the inner wall as seen on the left part of Fig. [3] The 
structure implodes and radial forces shift particles closer to the inner wall. In a final stage, only 
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Integrated density N(rft) 



Integrated density N(rft) 



M 





8 10 



8 10 



Fig. 3. Time evolution of the azimuthally integrated charge density N in the plasma column, on the 
left for m - 3 and on the right for m - 7. At the end of the runs, particles stay very close to the 
inner wall. 



few particles have moved outwards while the main fraction drifted close to the inner boundary. 
Therefore we do not observe any significant diffusion radially outwards. 

The situation is very similar for the second run we show, namely A2. Here also, the linear 
regime exists until f a; 60 and the observed growth rate is in close agreement with the linear 
analysis (red straight line on the right part of Fig.[Tll. The plasma column stays roughly confined 
between Ri and R2, right part of Fig. [3] no strong density rearrangement is perceptible during 
this linear phase. During the transition from the linear to the non-linear state, the m - 1 structure 
emerges in the density map as shown on the right part of Fig. ID Vacuum gaps are formed and seven 
vortices rotate around the cylinder axis until some of them merge together, destabilizing the whole 
system and shrinking down to the inner wall. 

These conclusions apply to all of our runs starting with an initial equilibrium as a plasma an- 
nulus. Such behavior is expected because of the confinement theorem which states that on average, 
the radial excursion of particles is lim ited due to angu lar momentum conservation of the system 
(particles + electromagnetic field), see iDavidsonI d 1 9901) . 



3.2. Electrosphere 

Next, we switch to the results obtained for the electrosphere. As in our previous works, the rotation 
profile is chosen to mimic the rotation curve existing in the 3D electrosphere. We remind that 
different analytical expressions for the radial dependence of £1 are chosen by mainly varying the 
gradient in differential shear as follows 



Q(r) = (2 + tanh[a' (r - ro)] e'^' ) 



(26) 



Q» is the neutron star spin and r is normalized to the neutron star radius, R,. The values used are 
Usted in Table |2] In all cases, Q. starts from corotation with the star Q = O, at r = 1, followed by a 
sharp increase around r - 6 for £li 2 and a less pronounced gradient around r = 10 for ^3. Finally 
the rotati on rate asympt otes twice the neutron star rotation speed for large radii, see discussion and 
figures in iPetril (l2007bl) . In this paper, we show results for the rotation profiles Qi and Q2 and 
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Table 2. Parameters for the three rotation profiles used to mimic the azimuthal velocity of the 
plasma in the electrospheric disk. 



n 


a 


f} 


'0 




3.0 


5 X 10-5 


6.0 




1.0 


5 X 10-5 


6.0 




0.3 


5 X 10-5 


10.0 



Model Mode m 





3 


6.669 + 1.8715 


i 


ill 


8 


16.74 + 3.1925 


i 


ill 


15 


30.98 + 1.6313 


i 


^2 


3 


6.761 + 0.9759 


i 


£22 


5 


10.90 + 0.765 


i 




3 


5.643 + 0.3394 


i 



Table 3. Eigenvalues co for the electrospheric model Q.i, Q,2 and Q3 for some azimuthal modes m. 
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Fig. 4. Evolution of the total electrostatic increase in energy for the differentially rotating plasma 
column obtained from our 2D PIC simulations (black dots) and compared with the linear growth 
rate obtained from Tab.[3](red line), on the left for the case Q.i,m - 3 and on the right for Q.2, m - 3. 
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Fig. 5. Snapshot of the charge density corresponding to the runs shown in Fig.|4]at the transition 
phase between linear and non-linear regime. The m = 3 patterns are clearly identified for suffi- 
ciently small times. 
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Fig. 6. Time evolution of the azimuthally integrated charge density TV associated with the runs 
shown in Fig.|4] 

a perturbation pattern m - 3. Note that the plasma is entirely filling the cylinder from Wi to W2 
without vacuum gaps. First, the evolution of the electrostatic energy is shown in Fig. |4] In both 
cases, the observed growth rates during the linear stage are recovered, in agreement with our linear 
analysis (compare the black dots with the red straight line) for f < 8 in the first run and f < 12 
in the second case. After this first stage, a transition to non-linear regime appears. The density 
perturbation becomes significant and is clearly recognizable on the density map. Fig. |5] Because 
of small irregularities, the 3 vortices have different rotational speeds, some of them overtake the 
slowest one and coalesce to a larger unique vortex. Eventually, the last vortex is overtaken and in 
the final situation, only one vortex survives and drags towards the neutron star surface (the inner 
wall). Indeed, inspecting Fig. |6] the highest charge densities are found close to 7J 1 . As a conclusion, 
here again, we found no significant charge escaping the system and trying to reach the outer wall. 

To summarize these simulations, we demonstrated that our 2D PIC code is able to accurately 
reproduce the growth rates of the diocotron instability in the linear stage of its development in 
accordance with our previous analytical study. 

The long time evolution in the non-linear stage showed no significant particle diffusion across 
magnetic field lines. Current flowing in the system is not efficient. 

Nevertheless, the pulsar magnetosphere is believed to be subject to copious pair creation feed- 
ing the electrosphere with highly relativistic freshly born electrons and positrons. This external 
source of charge can drastically change the long time behavior of the diocotron instability. That 
is why in the next and last paragraph, we add some particle injection mechanism into the plasma 
column (or electrosphere) in order to represent in a simple manner the pair cascade phenomenon. 

4. RESULTS WITH PARTICLE INJECTION 

The most interesting case is certainly the one assuming that the non-neutral plasma is not iso- 
lated but an open system from which particles can freely enter or exit. Our main focus is on the 
electrospheric plasma around a neutron star Close to its surface, due to the high magnetic field in- 
tensity, some quantum electrodynamical processes can convert photons into electron-positron pairs. 
Therefore, in a crude way, we can think of it as a (time-dependent) source feeding the electrosphere 
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Fig. 7. Time evolution of the azimuthally integrated charge density N in the plasma column. On 
the left, a slow increase in the particle radial extension is visible, whereas for the case on the right, 
a longer time scale is need to clearly recognize this effect. 

with charged particles. Moreover, these charges are responsible for synchrotron and curvature ra- 
diation. All these phenomena can drastically affect the behavior of the diocotron instability. 

Our goal in this last section is to study in detail the first mentioned effect, namely an external 
source of charges feeding the plasma with a flux of particles. In the most general picture, the spatial 
and time behavior can be quite complicated. To simplify, we will assume that the source is located 
at an arbitrary fixed radius and that the particle flux is constant in time. 

In order to investigate the consequences of particle injection into the system, we run exactly the 
same configurations as above but add a source of particle which are injected at a radius chosen to 
be r = 2 Wi. Some simulation results are described in the following subsections. 

4. 1 . Plasma column 

First, we consider again the plasma annulus with constant charge density and trapped between the 
two conducting walls. In all these runs, we adopt exactly the same conditions as before except that 
we remove the initial density perturbation. Indeed, we are not interested in the linear growth rate 
anymore. In this way, no matter what is the exact pattern of the perturbation. 

The integrated density is shown in Fig. |2]and should be compared with Fig. [3] Note that the 
color scales are different in both situations. At the beginning of the runs, the plasma annulus rotates 
around its axis without any significant perturbation. However, as time goes, around f = 60 and 
f = 100 for the first and second run respectively, the equilibrium configuration has been destroyed 
and particles shrink down to the inner wall. As the simulations start, the location of the source 
is easily determined on these maps. Particles are injected at a radius r = 2Wi, depicted on the 
plots by a thin strip. Typical perturbations of the density within the annulus are shown in Fig. [8] 
In both cases, it seems that the fastest growing mode corresponds to a pattern m - 6. The source 
of particles corresponds to the inner annulus. In the first simulation, on the left figure, there exists 
a strong interaction between both rings because they show exactly the same azimuthal pattern, at 
least, at the beginning of the non-linear stage. However, this effect is much less pronounced in 
the second run, on the right, due to the larger separation between both these rings. From a careful 
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Fig. 8. Snapshot of the charge density in the plasma column having the same initial conditions as in 
Fig.|2] Due to particle injection, a second inner ring appears, with size slowly growing in time. The 
weak initial particle noise is sufficient to ignite the unstable modes which interact with the inner 
annulus. 

comparison between the cases with and without injection, we conclude that the radial extension at 
the end of the simulation is significantly larger (on averaged) when a source is present. For instance, 
for the annulus A 1 , the integrated density does not extend farther than r = 2 Wi when the source is 
switch off whereas its extension reaches r - 5W\ when charges are injected. The diffusion process 
would have been even more drastic if the time span would have been longer We will make this 
statement more clear in the two last sets of run. The annulus configuration was just a starting point 
to relate our numerical results to analytically known solutions for the growth rates. 

4.2. Electrosphere 

Next, we move on to the behavior of the electrosphere which is actually the main purpose of this 
work. We take the same configurations as in the previous section and add the source located at 
r - 2Wi. However, in order to look for particle diffusion farther away than the outer wall located 
at r = 20 Wi, we used an initial box larger with size W\ = 1 and W2 - 100. This allows us 
to investigate large radial excursion of the particle and strong indication that charges effectively 
diffuse across the magnetic field lines. 

Let us have a look on the charge density at the end of the runs. The final snapshot is given in 
Fig- m for ^^1 on the left, and Q.2 on the right. We recall that at the first place, the electrospheric 
plasma was only filled between r = Wi and r = 20 Wi . In the above pictures, we clearly see a 
non negligible charge density up to radii larger than (40 - 50) Wi or so. This is an unambiguous 
evidence of the efficiency of the diocotron instability to generate a current across the field lines. 
This fact is supported by inspecting Fig. [TOl The integrated density shows a monotonic and regular 
expansion. During a first rearrangement state, the plasma shrink down and approaches the neutron 
star surface. This is the meaning of the thin horizontal strip in green-yellow between r - Wi and 
r = 20 Wi for time less than roughly f < 20. Thus a gap is forming, but after sufficient particles 
have been injected into the system, these vacuum regions will be replenished and expand even 
farther radially outwards. 
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Fig. 9. Snapshot of the charge density for the differentially rotating electrosphere, on the left for 
the case Qi and on the right for Q2. 
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Fig. 10. Time evolution of the azimuthally integrated charge density N in the electrosphere for the 
case Qi on the left and for Q2 on the right. 

4.3. Vacuum 

Finally, it is also possible to start with the extreme situation containing no particle, i.e. vacuum. 
This is certainly the most relevant case to show the efficiency of particle transport across magnetic 
field lines. 

As in the previous cases, we inject particles at r = 2 Wi and let the system evolve in a self- 
consistent manner according to the electric drift approximation. In the first run, the two conducting 
walls are relatively close to each other, with Wi = 1 and W2 -20W\. In this way, we are able to see 
what is happening in the region close to the injection radius. Depending on the particle injection 
flux, the plasma column is growing radially outwards more or less quickly. As seen in the density 
snapshot shown on the left part of Fig. (TT] at the end of the simulation, particles fill almost all 
the space between the two walls. Charge transport across magnetic field lines seems to be very 
efficient. This is confirmed by inspecting the left part of Fig. [12] in which we see and monotonic 
increase in outer edge of the integrated charge density. Starting with r around 2 W\, the plasma 
annulus has non negligible density at r = 13 Wi at the end of the run. 
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Finally, in order to clearly pick out the long time behavior of this diocotron instability, we per- 
formed a last run with higher resolution Nr x - 400 x 512 as well as larger distance between 
the conducting walls, Wi = 1 and W2 - 100. Therefore we are in a configuration suitable to inves- 
tigate the non-linear long time evolution of the instability. The most relevant results are exposed on 
the right part of Fig. [TT] and Fig. [121 for the density and the integrated density respectively. Here 
again, for a sufficient long time, plasma fills the entire vacuum and spreads out in whole space. It 
is clearly seen that the barycenter of the charges drifts to larger and larger radii. This last run is 
an unambiguous proof demonstrating the crucial role of this non-neutral plasma instability for the 
modeling of pulsar magnetosphere. Plasma injection in the electrosphere will unavoidably push 
particles farther away from the neutron star. 




Fig. 11. Snapshot of the charge density at the final time of the simulation, starting with vacuum. 
4.4. Discussion 

We showed that the non-neutral plasma behaves very differently when isolated compared to the 
situation in which a source of particles is added. In the former case, charges tend to migrate ra- 




Fig. 12. Time evolution of the azimuthally integrated charge density N, starting with vacuum. 
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dially inwards whereas in the latter case, the total charge of the system is not conserved and they 
are allowed to move radially outwards. This migration to larger distances is induced by the in- 
creasing electric repulsion between particles. In the initial (unstable) equilibrium state, the Lorentz 
force, q E + V A B, compensates exactly the centrifugal force. The self-field parameter, defined 
by Se = luj^/Q.^ - Ingm^lB^, is such that an equilibrium is permitted. For a constant den- 
sity plasma column it corresponds to < 1 . We introduced the plasma frequency (squared) by 
(jjp = «e e^/we eo, the cyclotron frequency by Q^. = e B/nif., the mass of the particle mg (electron or 
positron in our case) and the particle density number «£■ 

Finally, due to the evident dichotomy between an isolated system for which migration inwards 
is observed and a system with injection of particles and a radially outwards expansion of the plasma 
as consequences, we could imagine a situation in between such that no significant radial motion 
is induced. However, we do not believe that such a regime is reachable, because a charge density 
increase will inevitably lead to a growing self-field parameter Sg to values eventually larger than 
those expected to find an equilibrium solution (which is unity for the constant density plasma 
column). Remember that the self-field parameter is by definition proportional to the particle density 
number for a fixed magnetic field intensity. For a plasma column, it is well known that the limit Se = 
1 corr esponds to the Brillouin zone and forbids radial confinement above this limit. 



Davidson 



(Il990l) . For this reason, there will be no particle injection rate at which both effects annihilate 
exactly. The secular evolution is always particle diffusion across magnetic field lines to larger 
distances. 

A qualitative and illustrative picture is as follows. For simplicity we take the example of the 
plasma column without inner voids, uniform density and a constant uniform axial magnetic field. 
The rotation profile for this column is therefore constant, in other words, the plasma column is 
in solid body rotation. Because the plasma is confined, the self-field parameter should be less 
than unity. If some particles were added to the system, keeping the magnetic field constant, the 
particle number density «e increases and so the self-field parameter with the same proportion. 
When the Brillouin limit is reached, confinement is broken and particles diffuse to larger radii. A 
new equilibrium state could be reached in principle whenever becomes less than unity due to 
particle dilution in a greater volume. Because the timescale for charge diffusion is much less than 
the azimuthal rotation period, we can think about a slow almost adiabatic motion. The column 
always readjusts to a new "quasi"-equilibrium state. Thus, electrostatic repulsion is at the heart of 
the diffusion process. 

Finally, we give an estimate for the time evolution of the outer boundary of the system in the 
following way. In our simulations, at a fixed injection rate, the radial excursion slows down because 
particles have more space to "dilute" in (area in 2D (or volume in 3D) proportional to 2 n rdr{ dz)) 
and therefore decreasing more efficiently their density number as r grows. This is indeed observed 
in our runs. The outer edge of the plasma filled region /?edge can roughly be fitted as evolving in 
time as a function of ^ft. Indeed, assuming that the number of particles injected per unit time is T, 
they should be contained within the radius /?edge- Assuming an uniform density within the column, 
starting from r = at f = 0, we get T t (dz) ~ rtgn R^, ( dz) leading to the expected yfi function 



(27) 
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This curve is shown as a black solid parabola line on Fig. [10] and [12] The constant factor in front of 
the ^ft was adjusted individually for each run. 



5. CONCLUSION 

We designed a 2D electrostatic PIC code to study the long time non-linear evolution of the dio- 
cotron instability in a background electromagnetic field with a particular emphasize on the pulsar 
magnetosphere where we believe that non-neutral plasmas made of pairs form an electrosphere 
known to be diocotron unstable. 

We first check our PIC code by looking for the linear growing stage of the diocotron unstable 
modes. We found very good agreement between the analytical predictions from the linear analysis 
and the growth rates derived from our PIC code. 

These large growth rates quickly destroy the laminar flow pattern due to differential rotation. 
Fluid elements merge together to form some macroparticles with high density and leaving large 
vacuum space between them. Moreover, if the system remains isolated, there is no significant par- 
ticle transport across magnetic field lines. Nevertheless, we found that when adding some external 
source of charge, there is a net outflow of charge flowing towards the light cylinder. This source is 
believed to come from pair creation in the innermost part of the electrosphere. Non-neutral plasma 
instabilities like the one presented here are an efficient way to transport particles from regions close 
to the neutron star to the region where the wind is launched. The connection between the wind and 
the close magnetosphere still needs to be explain. We think that our work will be a fruitful alterna- 
tive to existing models. 

In this approach, we neglected the magnetic field perturbation induced by the plasma. This 
assumption is true only for plasma flowing far away from the light cylinder (and inside it). In order 
to take into account the magnetic field back reaction caused by the plasma, we need to solve the full 
set of Maxwell equations . This work w ill be presented in a forthcoming paper and will extend the 
results already founded in lPetril (l2007ah about some relativistic stabilization effects on the diocotron 
instability. 
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